Passive Tracking of Underwater Acoustic Sources with Sparse Innovations

ABSTRACT

A system and method involve acoustic source localization using passive sonar and capitalizing on the sparse nature of a source location map (SLM). Two types of sparsity are exploited, namely sparsity in the support of the SLMs and sparsity in the innovations across consecutive SLMs. The first type is motivated by the desire to construct SLMs whose non-zero entries corresponded to locations where sources are present. The second type of sparsity is motivated by the observation that few changes occur in the support of consecutive SLMs. Per time instant, an SLM is obtained as the solution of a regularized least-squares problem, where the regularization terms are chosen to encourage the desired sparse structures in each SLM and the innovations. Each SLM may be obtained via a specifically-tailored, computationally-efficient proximal gradient algorithm.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a continuation-in-part of commonly-assigned U.S.patent application Ser. No. 14/285,400, entitled “Multitask LearningMethod for Broadband Source-Location Mapping of Acoustic Sources”, thecontent of which is fully incorporated by reference herein.

FEDERALLY-SPONSORED RESEARCH AND DEVELOPMENT

This invention is assigned to the United States Government and isavailable for licensing for commercial purposes. Licensing and technicalinquiries may be directed to the Office of Research and TechnicalApplications, Space and Naval Warfare Systems Center, Pacific, Code72120, San Diego, Calif., 92152; voice (619) 553-5118; emailssc_pac_T2@navy.mil; reference Navy Case Number 103809.

BACKGROUND

Tracking acoustic sources via passive sonar is a challenging task commonto several underwater monitoring and surveillance systems. Classicaltracking approaches based on matched-field tracking and Kalman filteringtechniques are impractical due to their large computational and storagerequirements.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a block diagram of an embodiment of a system in accordancewith the Passive Tracking of Acoustic Sources with Sparse Innovations.

FIG. 2 shows a diagram of an embodiment of a data collection model usedin accordance with the Passive Tracking of Underwater Acoustic Sourceswith Sparse Innovations.

FIG. 3 shows a flowchart of an embodiment of a replica generation modelin accordance with the Passive Tracking of Acoustic Sources with SparseInnovations.

FIG. 4 shows a diagram illustrating an example of an operationalenvironment and corresponding source location map.

FIG. 5 shows a diagram of a broadband source localization map.

FIG. 6 shows a diagram of the structure of the regression coefficientmatrix S.

FIG. 7 shows a diagram illustrating the structure of the real-valuedmatrix {hacek over (S)}, which comprises the real and imaginary parts ofthe complex-valued matrix S as two different blocks.

FIG. 8 shows a graph illustrating the concept of a majorizer used by theproximal gradient method.

FIG. 9 shows a diagram of an ADMM algorithm for computing the proximalgradient update for pairs of rows in {hacek over (S)} that take nonzerovalues.

FIG. 10 shows a diagram of a PG algorithm developed for estimatingS_(m).

FIG. 11 shows a diagram of an accelerated PG algorithm developed forestimating S_(m).

FIG. 12 shows a flowchart of an embodiment of a method in accordancewith the Passive Tracking of Acoustic Sources with Sparse Innovations.

FIG. 13 shows a diagram of the propagation model used to characterizethe propagation environment of SWellEX-3.

FIG. 14 shows the source range (bottom panel) and the water depth (toppanel) along a trajectory of a towed acoustic source during SWellEX-3.

FIG. 15 shows a diagram illustrating an estimated depth track obtainedfor the source trajectory shown in FIG. 14 using matched field tracking.

FIG. 16 shows a diagram illustrating an estimated depth track obtainedfor the source trajectory shown in FIG. 14 using an embodiment of amethod in accordance with the Passive Tracking of Acoustic Sources withSparse Innovations.

FIG. 17 shows a diagram illustrating an estimated range track obtainedfor the source trajectory shown in FIG. 14 using matched field tracking.

FIG. 18 shows a diagram illustrating an estimated range track obtainedfor the source trajectory shown in FIG. 14 using an embodiment of amethod in accordance with the Passive Tracking of Acoustic Sources withSparse Innovations.

FIG. 19 shows a diagram illustrating a Source Location Map (SLM)obtained using an accelerated proximal gradient (PG) solver.

FIG. 20 shows a graph illustrating faster convergence using anaccelerated PG solver versus a PG solver.

FIG. 21 shows a diagram illustrating range tracks obtained for twosources using an embodiment of a method in accordance with the PassiveTracking of Acoustic Sources with Sparse Innovations.

DETAILED DESCRIPTION OF SOME EMBODIMENTS

Reference in the specification to “one embodiment” or to “an embodiment”means that a particular element, feature, structure, or characteristicdescribed in connection with the embodiments is included in at least oneembodiment. The appearances of the phrases “in one embodiment”, “in someembodiments”, and “in other embodiments” in various places in thespecification are not necessarily all referring to the same embodimentor the same set of embodiments.

Some embodiments may be described using the expression “coupled” and“connected” along with their derivatives. For example, some embodimentsmay be described using the term “coupled” to indicate that two or moreelements are in direct physical or electrical contact. The term“coupled,” however, may also mean that two or more elements are not indirect contact with each other, but yet still co-operate or interactwith each other. The embodiments are not limited in this context.

As used herein, the terms “comprises,” “comprising,” “includes,”“including,” “has,” “having” or any other variation thereof, areintended to cover a non-exclusive inclusion. For example, a process,method, article, or apparatus that comprises a list of elements is notnecessarily limited to only those elements but may include otherelements not expressly listed or inherent to such process, method,article, or apparatus. Further, unless expressly stated to the contrary,“or” refers to an inclusive or and not to an exclusive or.

Additionally, use of the “a” or “an” are employed to describe elementsand components of the embodiments herein. This is done merely forconvenience and to give a general sense of the invention. This detaileddescription should be read to include one or at least one and thesingular also includes the plural unless it is obviously meantotherwise.

The embodiments disclosed herein use sparse-signal modeling tools todevelop a broadband tracking algorithm for multiple underwater acousticsources that effectively tracks the broadband source-location map (SLM)as it evolves over time. Spectral data from multiple frequency bands areprocessed ‘coherently’ so as to unambiguously agree on the sourcelocations across frequencies. The tracking problem is cast so that thesparsity inherent in the SLM and in the SLM-innovations can beexploited. A numerical solver based on the proximal gradient (PG) methodand the alternating directions method of multipliers (ADMM) is developedfor SLM estimation.

Passive sonar enables monitoring and surveillance systems to operatewithout radiating sound into the water; hence, it is often employed inapplications where concealment and low environmental impact are desired.Acoustic data collected over time can be used for sketching sourcetracks by, for example, plotting source-location estimates over time.Tracking capitalizes on the temporal structure inherent to sourcetracks, which are always constrained by the kinematic features of thesource, to improve source-location estimates. However, using classicaltracking methods, such as Kalman filtering or Matched-field tracking(MFT), to develop a passive acoustic tracker poses significantcomputational challenges.

MFT is a passive acoustic tracking approach that builds on a generalizedunderwater beamforming technique called matched-field processing (MFP).MFP postulates a grid of tentative source locations and relies on anacoustic propagation model to obtain model-predicted acoustic pressures,hereafter referred to as replicas, at an array of hydrophones. Replicasare “matched” to acoustic measurements collected by an array ofhydrophones to construct a surface that summarizes the acoustic-powerestimates across each of the grid locations.

MFT relies on a sequence of ambiguity surfaces obtained at consecutivetime intervals for constructing tracks. Note that these ambiguitysurfaces are constructed independently from each other withoutexploiting their temporal correlation. Then, MFT constructs a graphconnecting grid points on consecutive ambiguity surfaces. Each possiblepath over this graph connecting the initial and final ambiguity surfaceis scored based on, e.g., its average ambiguity surface value, and thepath with the largest score serves as a track estimate. Unfortunately,the complexity of MFT grows quickly with the size of the grid and thenumber of ambiguity surfaces used. Constraints obtained from thesource's kinematics are often used to limit the number of tracks to beexplored.

Sparsity-driven Kalman-filter approaches can be used for trackingacoustic sources over a grid. In these approaches, the entire grid takesthe place of the state variable. It is presumed that only a few entriesof the grid, that is, those corresponding to the locations of thesources, take nonzero values. Thus, the state variable is considered tobe sparse. Unfortunately, the high dimensionality of the grid rendersimpractical any tracker that computes the full covariance matrix of thestate variable.

The broadband tracking approach disclosed herein builds on thesparsity-driven framework outlined above. It aims to construct sparsesource location maps (SLMs), one per frequency, while exploiting theirtemporal dependence. Only those grid points whose locations correspondto a source location take a nonzero value, which represents thecomplex-valued acoustic signature of the source. The tracker guaranteesthe support of the various grids to coincide, thereby assuringunambiguous source location estimates across frequencies. Different fromKalman-based trackers, our tracker relies on prior SLM estimates only.

An innovation is defined as a change in the support of consecutive SLMs.The proposed tracker controls the number of innovations allowed tooccur. It is assumed that innovations are sparse, and thus few sourceschange their location between consecutive SLMs. The disclosedembodiments involve a broadband tracker in which source-locationinformation is shared across frequencies. The resulting complex-valuedoptimization problem features a compounded regularizer that encouragessparse SLMs and sparse innovations. Iterative solvers based on acombination of the PG method and ADMM are developed.

FIG. 1 shows a block diagram of an embodiment of the operational conceptof a system in accordance with the Passive Tracking of Acoustic Sourceswith Sparse Innovations. System 10 includes K acoustic sources 20, 30,and 40 each radiating sound underwater. Although the acoustic sourcesare presumed to be mobile, thus justifying the dependence of theirlocations {r_(k)(τ)}_(k=1) ^(K) on the time τ ∈

, no assumptions about their kinematics are made. Each r_(k)(τ) ∈

^(d) is given in cylindrical coordinates comprising the source's range,depth (with respect to the sea surface), and azimuth, with d ∈ {1,2,3}.

An acoustic sensor array 50 detects acoustic signals 22, 32, and 42 fromsources 20, 30, and 40, respectively. As an example, an acoustic sensorarray 50 is an array of N hydrophones with known and arbitrary geometry,where N>0. Array 50 is used to collect a time series of acousticpressure vectors {y (τ) ∈

^(N): τ ∈

} with entries [y (τ)]_(n) ∈

denoting the acoustic pressure measured by the n-th hydrophone in thearray at time τ. Note that although the framework is agnostic to thespecific geometry of the array, its geometry affects the definition ofsource location. For instance, data gathered with an array that featureshorizontal aperture but no vertical aperture primarily providesinformation about the sources' azimuth only (d=1), although source rangeand depth can be determined depending on the closeness of the source tothe end fire direction, the length of the array, and the bandwidthconsidered. On the other hand, data collected on a vertical array withno horizontal aperture provide information on the source range and depthonly (d=2), unless the bathymetry is not flat. Data gathered with anarray featuring vertical and horizontal aperture provides informationabout the sources' range, depth, and azimuth (d=3).

Localization and tracking algorithms that use {y(τ): τ ∈

} directly are often challenged by the high sampling and computationalrequirements necessary to reconstruct the channel impulse responsesbetween the source locations and the hydrophones. Instead, this workdevelops a frequency-based passive-tracking approach that does notrequire estimating the channel impulse responses. To this end, {y(τ): τ∈

} is sequentially processed per hydrophone by computing short-timeFourier transforms (STFT) of length T₀. Consecutive blocks of theacoustic-pressure time series data can overlap and be scaled (perhydrophone) by a carefully chosen window function so as to reducedsidelobes and decrease the variability in the spectrum of theacoustic-pressure series due to noise. Fourier coefficients at Ffrequencies {ω_(f)}_(f=1) ^(F) across the N hydrophones are gathered inY^(m):=[y₁ ^(m), . . . , y_(F) ^(m)] ∈

^(N×F), where [Y^(m)]_(n,f) ∈

denotes the Fourier coefficient corresponding to ω_(f) for the n-thhydrophone in the m-th STFT block.

Each y_(f) ^(m) is modeled as

y _(f) ^(m)=Σ_(k=1) ^(K) s _(k,f) ^(m) p _(k,f)+ε_(f) ^(m) , ∀m, f,  (Eq. 1)

where s_(k,f) ^(m) ∈

denotes the Fourier coefficient at ω_(f) of the acoustic signature ofthe k-th source obtained by the m-th STFT, p _(k,f) ∈

^(N) the model-predicted Fourier coefficient vector at the array for thek-th source at frequency ω_(f) normalized so that ∥ p _(k,f)∥₂=1, andε_(f) ^(m) the spectral components of the array's measurement noise atω_(f) for the m-th STFT block. The replicas p _(k,f)'s are obtainedusing a model that characterizes the acoustic propagation environmentand the geometry of the array.

Given K, the goal of the spectral passive-acoustic tracking problem isto recursively estimate the locations {r_(k)(τ_(M))}_(k=1) ^(K) of theacoustic sources 20, 30, and 40 at τ_(M)=└(1−α)T₀M┘ based on thesequence of Fourier coefficient matrices {Y^(m)}_(m=0) ^(M), where α ∈[0,1) denotes percentage of overlap between consecutive STFT blocks.Note that here τ_(M) denotes the time index corresponding to thebeginning of the M-th temporal acoustic-data block. Even if all s_(k,f)^(m)'s were known, finding {r_(k)(τ_(M))}_(k=1) ^(K) based on (Eq. 1) ischallenging due to the non-linear relationship between each r_(k)(τ_(M))and its corresponding p _(k,f), a relationship that in most cases is notavailable in closed form.

In some embodiments, the estimation of the locations of the acousticsources is performed by a processor 52 embedded within acoustic sensorarray 50, with the results being transmitted, via signal 54, to receiver60. In some embodiments, acoustic sensor array sends, via signal 54, thesignal measurements to receiver 60, which performs the requiredprocessing.

As shown in diagram 100 of FIG. 2, acoustic sources 110 generateacoustic signals 120 that propagate through the propagation environmenttowards the acoustic sensor array 140, which is this case is ahydrophone array. Acoustic signals are collected by the hydrophone array140 as a time series 130 of acoustic pressures. The times series datacaptured by each hydrophone in the array 140 is transformed to theFourier domain via a short-time Fourier transform (STFT) processing step150 yielding Fourier coefficients. The Fourier coefficients obtainedacross hydrophones for a set of F frequencies {ω₁, . . . , ω_(F)} aregrouped to construct Fourier coefficient vectors y_(m,f) ∈

^(N) at the end of processing step 160.

Referring to FIG. 3, diagram 200 shows an embodiment of a model-basedprediction component, where a postulated source location 210 is input212 into a model-based predictor 220. Based upon an acoustic model, suchas an underwater acoustic model characterizing the true propagationenvironment, predictor 220 yields replica vectors 222. Replicas aremodel-based predicted vectors of Fourier coefficients that are used bythe Passive Tracking of Acoustic Sources with Sparse Innovations tomodel the y_(m,f)'s and determine the source localizations.

FIG. 4 shows a diagram 300 illustrating an example of an operationalenvironment and the desired SLM. As shown, an underwater environment 310includes a sensor array 320 disposed on the seafloor, a ship 330 on thesurface of the water, and submarines 340 and 350 located underwater. AnSLM 360 includes a marker 362 representing the location of sensor array320, a marker 364 representing the estimated location of ship 330, amarker 366 representing the estimated location of submarine 340, and amarker 368 representing the estimated location of submarine 350.

A model for the y_(f) ^(m)'s that alleviates the challenges associatedwith the nonlinearities inherent to (Eq. 1) is proposed next. Let

:={r_(g) ∈

^(d)}_(g=1) ^(G), with G>>max{N, KF}, denote a grid of tentative sourcelocations over the region of interest. Each y_(f) ^(m) is now modeled as

y _(f) ^(m)=Σ_(g=1) ^(G) s _(g,f) ^(m) p _(g,f)+ε_(m,f) , ∀f,   (Eq. 2)

where s_(g,f) ^(m) denotes the unknown Fourier coefficient associated tothe acoustic signature at frequency ω_(f) for a source located at r_(g),p_(g,f) ∈

^(N) is the replica for ω_(f) corresponding to a source located at r_(g)normalized so that ∥p_(g,f)∥₂=1 and ε_(m,f) denotes the Fouriercoefficients at ω_(f) corresponding to the noise in the m-th block.

It should be noted that (Eq. 2) tacitly assumes that the acousticsources are located exactly on some r_(g) ∈

. Since G>>KF, most of the s_(g,f) ^(m)'s are expected to be zero at anygiven m. Only coefficients s_(g,f) ^(m) that correspond to the locationof the acoustic sources take nonzero values, and thus theircorresponding replicas participate in (Eq. 2). All s_(g,f) ^(m)'sassociated to locations where acoustic sources are absent are set tozero.

From the vantage point of (Eq. 2), finding estimates for{r_(k)(τ_(M))}_(k=1) ^(K) is tantamount to identifying the locations in

corresponding to the nonzero s_(g,f)^(m‘'s. Let (•)′ denote the transpose operator and s) _(f)^(m):=[s_(1,f) ^(m), . . . , s_(G,f) ^(m)]′ ∀f. Once an estimate forS^(m):=[s₁ ^(m), . . . , s_(F) ^(m)] ∈

^(G×F) is available, one can construct the broadband SLM. This can bedone by, e.g., plotting the pairs (r_(g), ∥

_(g) ^(m)∥₂) for all r_(g) ∈

, where

_(g) ^(m):=[s_(g,1) ^(m), . . . , s_(g,F) ^(m)]′ ∈

^(F) comprises the entries of the g-th row of S_(m). Source locationestimates

indexed by the index-set

⊂ { 1, . . . , K} are given by the locations that correspond to theK-largest entries in the map, that is

∈    k  q ( Eq .  3 )

A localization algorithm that exploits the inherent sparse structure ofS^(m) was previously proposed by Forero. This approach estimates S^(m)using {y_(f) ^(m)}_(f=1) ^(F) only, while enforcing a common-supportacross all its columns {s_(f) ^(m)}_(f=1) ^(F). The requirement on thesupport of S_(m) is justifiable since it is assumed that the acousticsignal radiated by each source spans all {ω_(f)}_(f=1) ^(F).Nevertheless, such an estimator for S_(m) does not capture the temporaldependencies inherent to source motion that are caused by the physicalconstraints on the kinematics of the acoustic sources.

An iterative estimator for S^(m) is disclosed herein. The distinctivefeature of this estimator is that it uses the previously estimatedS^(m−1) to capture the temporal evolution of the source locations. Atτ_(m), S^(m) is obtained as the solution of the following regularizedleast-squares problem

 1 2  ∑ f = 1 F   y F m - P f  s f  2 2 + ∑ g = 1 G  [ μ   g 2 + λ   g - g m - 1  2 ] , ( Eq .  4 )

where S:=[s₁, . . . , s_(F)] ∈

^(G×F),

′_(g) is the g-th row of S, P_(f):=[p_(1,f), . . . , p_(G,f)] ∈

^(N×G) is the matrix of replicas for ω_(f), and μ, λ>0 are tuningparameters. The regularization term in (Eq. 4) encourages both groupsparsity on the rows of S^(m) and sparsity in the innovations, that is,changes in the support between the SLMs at τ_(m−1) and τ_(m).

The first term of the regularizer promotes the common support for thecolumns of S^(m). The tuning parameter μ controls the number of nonzerorows in S^(m). The second term of the regularizer captures temporalinformation related to the previous SLM via the

_(g) ^(m−1)'s. The number of innovations allowed to occur betweenτ_(m−1) and τ_(m) is controlled via λ.

FIG. 5 shows a diagram 400 of an embodiment of a broadband SLM 410comprising a plurality of tentative locations (grid points) and anestimated source location 420 for q=2. Grid point 420 corresponds to apair (r_(g), ∥

_(g) ^(m)∥₂) that subsumes the corresponding location r_(g) and the sizeof the corresponding regression coefficients

_(g) associated to that location across frequencies {ω_(f)}_(f=1) ^(F).Although only one estimated source location 420 is shown for SLM 410,other SLMs may contain additional estimated source locations 420.

FIG. 6 shows a diagram 500 of the structure of the regressioncoefficient matrix S 510. Once an estimate for S has been obtained, itscolumns can be used to construct SLMs over

per ω_(f). Furthermore, a broadband SLM can be constructed using wholerows

_(g) ∈

^(F) 520 of S for each r_(g) ∈

. For instance, an SLM can be constructed by plotting the pairs (r_(g),∥

_(g)∥_(q)) for all r_(g) ∈

and q>1.

As stated in the below proposition, (Eq. 4) can be written as areal-valued convex optimization problem after representing allcomplex-valued variables by the direct sum of their real and imaginaryparts. Before stating the proposition, the following notation isintroduced {hacek over (y)}_(f) ^(m):=[Re(y_(F) ^(m))′, Im(y_(F)^(m))′]′, {hacek over (s)}_(f) ^(m):=[Re(s_(F))′, Im(s_(F))′]′, {hacekover (S)}:=[{hacek over (s)}₁, . . . , {hacek over (s)}_(F)], and

$\begin{matrix}{{\overset{\Cup}{P}}_{f}:=\begin{bmatrix}{{Re}\left( P_{f} \right)} & {- {{Im}\left( P_{f} \right)}} \\{{Im}\left( P_{f} \right)} & {{Re}\left( P_{f} \right)}\end{bmatrix}} & \left( {{Eq}.\mspace{14mu} 5} \right)\end{matrix}$

where Re(•) (Im(•)) denotes the real-part (imaginary-part) operator.Matrix {hacek over (S)} can be alternatively viewed in terms of its rowsas {hacek over (S)}=[{hacek over (

)}₁, . . . , {hacek over (

)}_(2G)]′ where the first (last) G rows correspond to the real(imaginary) parts of the rows of S.

The minimum of (Eq. 3) is equal to that of the following convexoptimization problem

$\begin{matrix}{{{\arg \frac{1}{2}{\sum\limits_{f = 1}^{F}{{{\overset{\Cup}{y}}_{F}^{m} - {{\overset{\Cup}{P}}_{f}{\overset{\Cup}{s}}_{f}}}}_{2}^{2}}} + {\mu {\sum\limits_{g = 1}^{G}{v_{g}}_{2}}} + {\lambda {\sum\limits_{g = 1}^{G}{{{\overset{\Cup}{v}}_{g} - {\overset{\Cup}{v}}_{g}^{m - 1}}}_{2}}}},} & \left( {{Eq}.\mspace{14mu} 6} \right)\end{matrix}$

where {hacek over (v)}_(g):=[{hacek over (

)}′_(g), {hacek over (

)}′_(g+G)]′ ∈

_(2F) ({hacek over (v)}_(g) ^(m−1)) corresponds to the direct sum of thereal and imaginary parts of

_(g) (

_(g) ^(m−1)) . FIG. 7 shows a diagram 600 illustrating the structure of{hacek over (S)}. It comprises two blocks that correspond to the realpart 610 of matrix S and the imaginary part 620 of matrix S. Theminimizer S^(m) of (Eq. 4) is given in terms of {hacek over (S)}^(m) asS^(m)={hacek over (S)}_(1:G) ^(m)+j{hacek over (S)}_(G+1:2G) ^(m), wherej:=√{square root over (−1)} and {hacek over (S)}_(g1:g2) ^(m)=[{hacekover (

)}_(g1), . . . , {hacek over (

)}_(g2)]′, for all g₁≦g₂, g₁, g₂ ∈ {1, . . . , G}.

Although (Eq. 6) is a convex optimization problem that can be solvedvia, e.g., interior point methods, such solver entails highcomputational complexity due to the large dimensionality of {hacek over(S)} and fails to exploit the sparse structure of {hacek over (S)}.Discussed below is a PG solver for (Eq. 5) that exploits its sparsestructure.

The problem in (Eq. 6) can be written asmin_({hacek over (S)} h({hacek over (S)})+θ({hacek over (S)}), where)

${h\left( \overset{\Cup}{S} \right)}:={\frac{1}{2}{\sum\limits_{f = 1}^{F}{{{\overset{\Cup}{y}}_{F}^{m} - {{\overset{\Cup}{P}}_{f}{\overset{\Cup}{s}}_{f}}}}_{2}^{2}}}$

denotes the continuously differentiable portion of the cost, andθ({hacek over (S)}):=μΣ_(g=1) ^(G)∥{hacek over (v)}_(g)∥₂+λΣ_(g=1)^(G)∥{hacek over (v)}_(g)−{hacek over (v)}_(g) ^(m−1)∥₂ thenon-differentiable portion of the cost. Note that the gradient ofh({hacek over (S)}) is Lipschitz continuous with Lipschitz constantL_(h):=max_(f=1, . . . , F) σ_(max)(P′_(f)P_(f)), whereσ_(max)(P′_(f)P_(f)) denotes the largest singular value of P′_(f)P_(f).That is, ∥∇h({hacek over (S)}₁)−∇h({hacek over (S)}₂)∥₂≦L_(h)∥{hacekover (S)}₁−{hacek over (S)}₂∥_(F), where ∇h({hacek over (S)}_(l))denotes the gradient of h with respect to {hacek over (S)} evaluated at{hacek over (S)}_(l).

The PG method can be interpreted as a majorization-minimization methodrelying on a majorizer H({hacek over (S)}; Z) for h, where Z:=[z₁, . . ., z_(F)] ∈

^(2G×F) is an auxiliary matrix. The majorizer H satisfies: (i) H({hacekover (S)}; Z)≧h({hacek over (S)}), ∀({hacek over (S)}); and, (ii)H({hacek over (S)}; Z)=h({hacek over (S)}) for Z={hacek over (S)}. Thespecific H used is

$\begin{matrix}{{H\left( {\overset{\Cup}{S};Z} \right)}:={{h(Z)} + {\sum\limits_{f = 1}^{F}{{\nabla{h_{f}\left( z_{f} \right)}^{\prime}}\left( {{\overset{\Cup}{s}}_{f} - z_{f}} \right)}} + {\frac{L_{h}}{2}{{\overset{\Cup}{S} - Z}}_{F}^{2}}}} & \left( {{Eq}.\mspace{14mu} 7} \right)\end{matrix}$

where h_(f)({hacek over (s)}_(f)):=1/2∥{hacek over (y)}_(f) ^(m)−{hacekover (P)}_(f){hacek over (s)}_(f)∥₂ ², and ∇h_(f)(z_(f)) denotes thegradient of h_(f) with respect to {hacek over (s)}_(f) evaluated atz_(f). That (Eq. 7) satisfies conditions (i) follows from the fact thatthe gradient of h is Lipschitz continuous, and that it satisfies (ii)follows after setting Z={hacek over (S)} in (Eq. 7). FIG. 8 shows agraph 700 illustrating a majorizer H 710 for the smooth function h 720.Note that the line intersection 730 illustrates condition (ii) above,that is, H must touch at the point Z={hacek over (S)}. With j denotingthe iteration index, the PG algorithm iteratively solves

$\begin{matrix}{{{\overset{\Cup}{S}}^{m}\lbrack j\rbrack} = {\min\limits_{\overset{\Cup}{S}}{\left\lbrack {{H\left( {\overset{\Cup}{S};{{\overset{\Cup}{S}}^{m}\left\lbrack {j - 1} \right\rbrack}} \right)} + {\theta \left( \overset{\Cup}{S} \right)}} \right\rbrack.}}} & \left( {{Eq}.\mspace{14mu} 8} \right)\end{matrix}$

From an algorithmic point of view, it is convenient to write H as afunction of the {hacek over (v)}_(g)'s. After performing some algebraicmanipulations on H and dropping all terms independent of {hacek over(S)}, (Eq. 8) can be written as

$\begin{matrix}{{{{\overset{\Cup}{S}}^{m}\lbrack j\rbrack} = {\min\limits_{\overset{\Cup}{S}}\left\lbrack {{\sum\limits_{g = 1}^{G}\; {\frac{L_{h}}{2}{{{\overset{\Cup}{v}}_{g} - {w_{g}^{m}\left\lbrack {j - 1} \right\rbrack}}}_{2}^{2}}} + {\theta \left( \overset{\Cup}{S} \right)}} \right\rbrack}},{where}} & \left( {{Eq}.\mspace{14mu} 9} \right) \\{{w_{g}^{m}\left\lbrack {j - 1} \right\rbrack}:={{{\overset{\Cup}{v}}_{g}^{m}\left\lbrack {j - 1} \right\rbrack} - {\left( \frac{1}{L_{h}} \right){d_{g}^{m}\left\lbrack {j - 1} \right\rbrack}}}} & \left( {{Eq}.\mspace{14mu} 10} \right)\end{matrix}$

is a gradient-descent step, with step-size

$\frac{1}{L_{h}},$

for the g-th row of {hacek over (S)}, and the entries of d_(g)^(m)[j−1], which correspond to those of the gradient of h_(f) withrespect to {hacek over (v)}_(g), are

$\begin{matrix}{\left\lbrack {D_{g}^{m}\left\lbrack {j - 1} \right\rbrack} \right\rbrack_{f} = \left\{ \begin{matrix}{{{- {\overset{\Cup}{p}}_{g,f}^{\prime}}{r_{f}^{m}\left\lbrack {j - 1} \right\rbrack}},} & {{f = 1},\ldots \mspace{14mu},F} \\{{{- {\overset{\Cup}{p}}_{{g + G},f}^{\prime}}{r_{f}^{m}\left\lbrack {j - 1} \right\rbrack}},} & {{f = {F + 1}},\ldots \mspace{14mu},{2\; F}}\end{matrix} \right.} & \left( {{Eq}.\mspace{14mu} 11} \right)\end{matrix}$

where r_(f) ^(m)[j−1]:={hacek over (y)}_(f)−{hacek over (P)}_(f){hacekover (s)}_(f) ^(m)[j−1]. Eq. 9 is often called the proximal operator ofθ with parameter

$\frac{1}{L_{h}}.$

Eq. 9 is decomposable across {hacek over (v)}_(g)'s. Per iteration j,the PG update in (Eq. 8) can be performed in parallel for every pair ofrows of {hacek over (S)} comprised in each {hacek over (v)}_(g) via

$\begin{matrix}{{{\overset{\Cup}{v}}_{g}^{m}\lbrack j\rbrack} = {\min\limits_{{\overset{\Cup}{v}}_{g}}{\left\lbrack {{\frac{L_{h}}{2}{{{\overset{\Cup}{v}}_{g} - {w_{g}^{m}\left\lbrack {j - 1} \right\rbrack}}}_{2}^{2}} + {\mu {{\overset{\Cup}{v}}_{g}}_{2}} + {\lambda {{{\overset{\Cup}{v}}_{g} - {\overset{\Cup}{v}}_{g}^{m - 1}}}_{2}}} \right\rbrack.}}} & \left( {{Eq}.\mspace{14mu} 12} \right)\end{matrix}$

The cost in (Eq. 12) is convex; however, it is non-differentiable due tothe compounded regularization term. This regularizer is such that aclosed-form update for {hacek over (v)}_(g) ^(m)[j] in general is notavailable. Thus, (Eq. 12) must be solved numerically while bearing inmind that the computational cost associated to each PG iteration hingeson that of solving (Eq. 12). Note that (Eq. 12) can be solved in closedform for the special case where {hacek over (v)}_(g) ^(m−1)=0_(2F),where 0_(2F) is a 2F×1 vector of zeros, and its solution in this case is

$\begin{matrix}{{{\overset{\Cup}{v}}_{g}^{m}\lbrack j\rbrack} = {{w_{g}^{m}\left\lbrack {j - 1} \right\rbrack}\left( {1 - \frac{\lambda + \mu}{L_{h}{{w_{g}^{m}\left\lbrack {j - 1} \right\rbrack}}_{2}}} \right)_{+}}} & \left( {{Eq}.\mspace{14mu} 13} \right)\end{matrix}$

where (•)₊=max{0,•}.

To gain further insight into the solution of (Eq. 12), when {hacek over(v)}_(g) ^(m−1)≠0_(2F), the Karush-Kuhn-Tucker (KKT) conditions may beused combined with the notion of subdifferential to characterize {hacekover (v)}_(g) ^(m)[j]. With V^(m)({hacek over (v)}_(g)) denoting thecost in (Eq. 12) and since (Eq. 12) is an unconstrained optimizationproblem, the KKT optimality conditions state that 0_(2F) ∈ ∂V^(m)({hacekover (v)}_(g) ^(m)[j]) where ∂V^(m)({hacek over (v)}_(g)) denotes thesubdifferential of V_(m) evaluated at {hacek over (v)}_(g). Thefollowing necessary and sufficient conditions for {hacek over (v)}_(g)^(m)[j] follow readily from the optimality condition.

When {hacek over (v)}g^(m−1)≠0_(2F), the following mutually-exclusiveconditions about (Eq. 12) hold

$\begin{matrix}{{{\overset{\Cup}{v}}_{g}^{m}\lbrack j\rbrack} = \left. 0_{2\; F}\Leftrightarrow{{{{{w_{g}^{m}\left\lbrack {j - 1} \right\rbrack} +} \propto_{g}^{m - 1}{\overset{\Cup}{v}}_{g}^{m - 1}}}_{2} \leq \frac{\lambda}{L_{h}}} \right.} & \left( {{{Eq}.\mspace{14mu} 14}a} \right) \\{{{\overset{\Cup}{v}}_{g}^{m}\lbrack j\rbrack} = \left. {\overset{\Cup}{v}}_{g}^{m - 1}\Leftrightarrow{{{{w_{g}^{m}\left\lbrack {j - 1} \right\rbrack} + {\beta_{g}^{m - 1}{\overset{\Cup}{v}}_{g}^{m - 1}}}}_{2} \leq \frac{\mu}{L_{h}}} \right.} & \left( {{{Eq}.\mspace{14mu} 14}b} \right)\end{matrix}$where ∝_(g) ^(m−1) =μ/L _(h) ∥{hacek over (v)} _(g) ^(m−1)∥₂ and β_(g)^(m−1)=1+λ/L _(h) ∥{hacek over (v)} _(g) ^(m−1)∥₂.

The above equations can be used to quickly screen whether {hacek over(v)}^(m)[j] equals 0_(2F) or {hacek over (v)}^(m−1). If neither of theseconditions is satisfied, then {hacek over (v)}^(m)[j] must be obtainedby solving (Eq. 12) numerically.

An efficient iterative solver for (Eq. 12) based on ADMM is developed inthis section. To this end, consider the following optimization problem,which is equivalent to (Eq. 12),

$\begin{matrix}{{{\min\limits_{{\overset{\Cup}{v}}_{g},\zeta_{g}}{\frac{L_{h}}{2}{{{\overset{\Cup}{v}}_{g} - {w_{g}^{m}\left\lbrack {j - 1} \right\rbrack}}}_{2}^{2}}} + {\mu {{\overset{\Cup}{v}}_{g}}_{2}} + {\lambda {\zeta_{g}}_{2}}}{{{{Subj}.\mspace{14mu} {to}}\mspace{14mu} \zeta_{g}} = {{\overset{\Cup}{v}}_{g} - {{\overset{\Cup}{v}}_{g}^{m - 1}.}}}} & \left( {{Eq}.\mspace{14mu} 15} \right)\end{matrix}$

The ADMM solver for (Eq. 15) relies on the augmented Lagrangian given by

$\begin{matrix}{{\mathcal{L}\left( {{\overset{\Cup}{v}}_{g},\zeta_{g},\gamma_{g}} \right)} = {{\frac{L_{h}}{2}{{{\overset{\Cup}{v}}_{g} - {w_{g}^{m}\left\lbrack {j - 1} \right\rbrack}}}_{2}^{2}} + {\mu {{\overset{\Cup}{v}}_{g}}_{2}} + {\lambda {\zeta_{g}}_{2}} + {\gamma_{g}^{\prime}\left( {\zeta_{g} - {\overset{\Cup}{v}}_{g} + {\overset{\Cup}{v}}_{g}^{m - 1}} \right)} + {\frac{n}{2}{{\zeta_{g} - {\overset{\Cup}{v}}_{g} + {\overset{\Cup}{v}}_{g}^{m - 1}}}_{2}^{2}}}} & \left( {{Eq}.\mspace{14mu} 16} \right)\end{matrix}$

where η>0 is a tuning parameter, and γ_(g) the Lagrange multipliervector associated to the equality constraint in (Eq. 15). With idenoting an iteration index, the ADMM iterations are

$\begin{matrix}{{{\overset{\Cup}{v}}_{g}\lbrack i\rbrack} = {\arg {\min\limits_{{\overset{\Cup}{v}}_{g}}{\mathcal{L}\left( {{\overset{\Cup}{v}}_{g},{\zeta_{g}\left\lbrack {i - 1} \right\rbrack},{\gamma_{g}\left\lbrack {i - 1} \right\rbrack}} \right)}}}} & \left( {{{Eq}.\mspace{14mu} 17}a} \right) \\{{\varsigma_{g}\lbrack i\rbrack} = {\arg {\min\limits_{\varsigma_{g}}{\mathcal{L}\left( {{{\overset{\Cup}{v}}_{g}\lbrack i\rbrack},\zeta_{g},{\gamma_{g}\left\lbrack {i - 1} \right\rbrack}} \right)}}}} & \left( {{{Eq}.\mspace{14mu} 17}b} \right) \\{{\gamma_{g}^{\prime}\lbrack i\rbrack} = {{\gamma_{g}^{\prime}\left\lbrack {i - 1} \right\rbrack} + {\eta\left( {{\zeta_{g}\lbrack i\rbrack} - {{\overset{\Cup}{v}}_{g}\lbrack i\rbrack} + {\overset{\Cup}{v}}_{g}^{m - 1}} \right.}}} & \left( {{{Eq}.\mspace{14mu} 17}c} \right)\end{matrix}$

The following equations show that updates to (Eq. 17a) and (Eq. 17b) canbe obtained in closed form via soft-thresholding, as they can beperformed in closed form as

$\begin{matrix}{{{\overset{\Cup}{v}}_{g}\lbrack i\rbrack} = {\frac{\rho_{g}^{m - 1}\left\lbrack {i - 1} \right\rbrack}{L_{h} + \eta}\left( {1 - \frac{\mu}{{{\rho_{g}^{m - 1}\left\lbrack {i - 1} \right\rbrack}}_{2}}} \right)_{+}}} & \left( {{{Eq}.\mspace{14mu} 18}a} \right) \\{{\zeta_{g}\lbrack i\rbrack} = {\frac{\chi_{g}^{m - 1}\left\lbrack {i - 1} \right\rbrack}{\eta}\left( {1 - \frac{\lambda}{{{\chi_{g}^{m - 1}\left\lbrack {i - 1} \right\rbrack}}_{2}}} \right)_{+}}} & \left( {{{Eq}.\mspace{14mu} 18}b} \right)\end{matrix}$where ρ_(g) ^(m−1) [i−1]:=η(ζ_(g) [i−1]+{hacek over (v)} _(g) ^(m−1) +L_(h) ·w _(g) ^(m) [j−1]+γ_(g) [i−1] and χ_(g) ^(m−1) [i−1]:=η({hacekover (v)} _(g) [i]−{hacek over (v)} _(g) ^(m−1))−γ_(g) [i−1].

The resulting ADMM algorithm is summarized as Algorithm 1 shown indiagram 800 of FIG. 9. Per iteration, its computational complexity isdominated by the evaluation of the Euclidean norms in (Eqs. 18a and18b), which require O(F) time complexity. With respect to theconvergence of Algorithm 1, note that the sequence {{hacek over(v)}_(g)[i]}_(i>0) does not need to converge to the optimal value of{hacek over (v)}_(g). Nevertheless, the results can be used to show thatevery limit point of {{hacek over (v)}_(g)[i]}_(i>0) corresponds to anoptimal solution of (Eq. 14), and that the sequence {∥ζ_(g)[i]−{hacekover (v)}_(g)[i]−{hacek over (v)}_(g) ^(m−1)∥₂}_(i>0) converges to zero,i.e., iterates {hacek over (v)}_(g)[i] and ζ_(g)[i] approach feasibilityas i→∞.

The resulting tracking PG algorithm is summarized as Algorithm 2 shownin diagram 900 of FIG. 10. In general, its per-iteration computationalcomplexity is dominated by the execution time of Algorithm 1 (line 12).Note, however, that Algorithm 1 is used only when {hacek over (v)}_(g)^(m)[j] ∉ {0_(2F), {hacek over (v)}_(g) ^(m−1)) and {hacek over (v)}_(g)^(m−1)≠0_(2F). Due to the high-sparsity regime in which Algorithm 2operates, few executions of Algorithm 1 are required. When Algorithm 1is not executed, its computational complexity is dominated by the updatein (Eq. 10) which entails O(NG) operations.

Algorithm 2 can be shown to converge to the solution of (Eq. 6) whilefeaturing a worst-case convergence rate of O(1/j). Thus its convergencemay be slow in practice, requiring up to several hundreds of iterationsto achieve a highly accurate solution. Recent works have shown that itis possible to improve the suboptimal convergence rate of the PG methodwhile maintaining its computational simplicity. These works propose todevelop accelerated PG algorithms that feature worst-case convergencerate of O(1/j²).

FIG. 11 shows a diagram 1000 of Algorithm 3, which is an accelerated PG(APG) algorithm developed for solving the problem in Eq. 6. Differentfrom Algorithm 2, the APG algorithm updates {hacek over (S)}^(m)[j]based on a slack variable Σ[j] constructed as a linear combination ofthe previous two estimates for {hacek over (S)}^(m). These update iscomputed similarly to that described for the PG algorithm as shown inFIG. 10. Note that the added time and space complexity of theaccelerated PG algorithm can be considered negligible, especially afterrecalling that both {hacek over (S)}^(m)[j] and {hacek over(S)}^(m)[j−1] are sparse.

FIG. 12 shows a flowchart of an embodiment of a method 1100 inaccordance with the Passive Tracking of Underwater Acoustic Sources withSparse Innovations. As an example, method 1100 may be performed bysystem 10, 100, and 200 as shown in FIGS. 1-3. Also, while FIG. 12 showsone embodiment of method 1100 to include steps 1110-1180, otherembodiments of method 1100 may contain fewer or more steps. Further,while in some embodiments the steps of method 1100 may be performed asshown in FIG. 12, in other embodiments the steps may be performed in adifferent order, or certain steps may occur simultaneously with one ormore other steps. Additionally, some or all of the steps of method 1100may be performed by processor 52 embedded within acoustic sensor array50, by receiver 60, or by other processing means operatively connectedto acoustic sensor array 50.

Method 1100 may begin with step 1110, which involves defining a grid

of G tentative locations r_(g), to localize more than one acousticsources, based on the location of an acoustic sensor array 50, where upto K acoustic sources 20, 30, and 40 are presumed to be located. In someembodiments step 1110 includes the step of estimating the number ofacoustic sources, while in other embodiments the number of acousticsources is predetermined.

Step 1120 involves using an acoustic model to compute, via themodel-based predictor 220, replicas that correspond to simulatedacoustic sources at locations in

, wherein the replicas are model-predicted Fourier coefficient vectorsat F frequencies {ω_(f)}_(f=1) ^(F) corresponding to the acousticpressure field 120 as sampled by an acoustic sensor array 140 having Nsensors. Step 1130 involves collecting, using acoustic sensor array 140,time series data 130 of actual acoustic measurements at each sensor ofthe acoustic sensor array caused by the acoustic sources 110.

Method 1100 may proceed to step 1140, which involves using a short-timeFourier transform (STFT) 150 on the collected time series data 130partitioned to m blocks to compute Fourier coefficient estimates atfrequencies {ω_(f)}_(f=1) ^(F) for all N sensors. Step 1150 involvesconstructing a set of Fourier coefficient vectors Y^(m):=[y₁ ^(m), . . ., y_(F) ^(m)] ∈

^(N×F), where [Y^(m)]_(n,f) ∈

denotes the Fourier coefficient corresponding to ω_(f) for the n-thsensor in the m-th STFT block, at step 160 using the Fouriercoefficients previously obtained via the STFT.

Step 1160 involves modeling Fourier coefficient vectors at ω_(f) for them-th measurement block of the collected time series data as y_(f)^(m)=Σ_(g=1) ^(G)s_(g,f) ^(m)p_(g,f)+ε_(m,f), ∀f, where s_(g,f) ^(m)denotes the unknown Fourier coefficient associated to the acousticsignature at frequency ω_(f) for a source located at r_(g), p_(g,f) ∈

^(N) is the replica for ω_(f) corresponding to a source located at r_(g)normalized so that ∥p_(g,f)∥₂=1, and ε_(m,f) denotes the Fouriercoefficients at ω_(f) corresponding to the noise in the m-th block.

Step 1170 involves setting to zero all values of s_(g,f) ^(m) associatedto tentative locations r_(g)where acoustic sources are presumed to beabsent. Step 1170 is an optional step that may or may not appear in allembodiments of the methods disclosed herein. It relies on so-calledpredictor screening rules that can identify with certainty points in

where acoustic sources are absent. A specific embodiment of apredictor-screening rule is described in the commonly-assigned U.S.patent application Ser. No. 14/285,400, entitled “Multitask LearningMethod for Broadband Source-Location Mapping of Acoustic Sources.”

Step 1180 involves obtaining, at time τ_(m), an estimate S_(m) as thesolution to

${{\min_{S \in {\mathbb{C}}^{G \times F}}{\frac{1}{2}{\sum\limits_{f = 1}^{F}\; {{y_{F}^{m} - {P_{f}s_{f}}}}_{2}^{2}}}} + {\sum\limits_{g = 1}^{G}\; \left\lbrack {{\mu {ϛ_{g}}_{2}} + {\lambda {{ϛ_{g} - ϛ_{g}^{m - 1}}}_{2}}} \right\rbrack}},$

where S:=[s₁, . . . , s_(F)] ∈

^(G×F),

′_(g) is the g-th row of S, P_(f):=[p_(1,f), . . . , p_(G,f)] ∈

^(N×G) is the matrix of replicas for ω_(f), and μ, λ>0 are tuningparameters. In some embodiments, the estimate of S^(m) in step 1180 isobtained via Algorithm 2 shown in FIG. 10. Step 1180 also involvesgenerating one or more SLMs over

per frequency ω_(f) using S^(m), wherein each location on a particularSLM is associated with the magnitude of its corresponding acoustic gainestimate |ŝ_(g,f)|, wherein estimates of the actual locations of the Kacoustic sources, {{circumflex over (r)}_(k)(τ_(m))}_(k∈κ), correspondto the locations of the K-largest coefficients |ŝ_(g,f)| depicted in theSLM. In some embodiments of step 1180, the SLMs are generated using aniterative problem solver based upon a PG method.

In some embodiments, step 1180 involves generating SLMs over

per frequency ω_(f) using each column s_(f) of S^(m) to construct theSLMs per frequency ω_(f) . In some embodiments, step 1180 involvesgenerating a broadband SLM over

comprising all frequencies ω_(f) used to compute S^(m) using a whole rowof S^(m) for each r_(g) ∈

. In some embodiments of step 1180, generating SLMs over

per frequency ω_(f) involves generating a broadband SLM over

using S^(m) by plotting the pairs (r_(g), ∥

_(g) ^(m)∥₂) for all r_(g) ∈

, where

_(g) ^(m):=[s_(g,1) ^(m), . . . , s_(g,F) ^(m)]′ ∈

^(F) comprises the entries of the g-th row of S^(m).

Method 1100 may be implemented as a series of modules, eitherfunctioning alone or in concert, with physical electronic and computerhardware devices. Method 1100 may be computer-implemented as a programproduct comprising a plurality of such modules, which may be displayedfor a user.

Various storage media, such as magnetic computer disks, optical disks,and electronic memories, as well as non-transitory computer-readablestorage media and computer program products, can be prepared that cancontain information that can direct a device, such as amicro-controller, to implement the above-described systems and/ormethods. Once an appropriate device has access to the information andprograms contained on the storage media, the storage media can providethe information and programs to the device, enabling the device toperform the above-described systems and/or methods.

For example, if a computer disk containing appropriate materials, suchas a source file, an object file, or an executable file, were providedto a computer, the computer could receive the information, appropriatelyconfigure itself and perform the functions of the various systems andmethods outlined in the diagrams and flowcharts above to implement thevarious functions. That is, the computer could receive various portionsof information from the disk relating to different elements of theabove-described systems and/or methods, implement the individual systemsand/or methods, and coordinate the functions of the individual systemsand/or methods.

The performance of the proposed broadband tracking algorithm isillustrated on the third Shallow-Water Evaluation Cell Experiment(SWellEX-3) dataset. The environment considered corresponds to the onein the third Shallow-Water Evaluation Cell Experiment (SWellEX-3), seediagram 1200 shown in FIG. 13. In SWellEX-3, a towed source transmittingat frequencies {53+16 k}_(k=0) ⁹ Hertz and a vertical line array 1210collecting acoustic data were used. In this analysis, only 9hydrophones, out of 64 hydrophones available, were used. Thesehydrophones were 11.25 m apart, having a total aperture of 90 m with thebottom element 6 m above the seafloor (water depth was 198 m).

A grid with G=20,000 locations spanning radial distances 0-10 km anddepths 0-198 m was used, as shown by the top portion of the diagram 1220ending at a depth of 198 m at line 1230. Depths between 198 m-228 m arerepresented by portion 1240, depths between 228 m-1028 m are representedby portion 1250, and depths below 1028 m are represented by portion1260. The grid's radial and vertical spacing were 50 m and 2 m,respectively. All replicas were computed with the KRAKEN normal-modepropagation model using the environmental model shown in FIG. 13.

Sample parameter values used in the model are: υ₁=1; 520 m/s, υ₂32 1;498 m/s, υ₃=1; 490 m/s, υ₄32 1;490 m/s, υ₅=1; 572 m/s, υ₆=1; 593 m/s,υ₇=1; 881 m/s, υ₈=3.246 m/s, υ₉5; 200 m/s, α_(b) ₁ =0.2 dB/m/kHz, α_(b)₂ =0.06 dB/m/kHz, α_(b) ₃ =0.02 dB/m/kHz, ρ_(b) ₁ =1.76 g/cm3, ρ_(b) ₂=2.06 g/cm3, and ρ_(b) ₃ =2.66 g/cm3.

Since ADMM is only used to update few of the v_(g)'s (in the order ofthe sparsity of the SLM) per PG iteration, its computational complexitydoes not significantly affect that of Algorithm 2. On the other hand,the PG algorithm takes a few hundreds of iterations to converge to anacceptable precision as defined by the quality of the source tracks.This motivates future work exploring both predictor screening rules andaccelerated PG methods as a mean to reduce the computational complexityof Algorithm 2 both by reducing G and the number of PG iterationsrequired.

The top portion the diagram 1300 shown in FIG. 14 illustrates thebathymetry 1310 along the range-trajectory of a towed acoustic sourceused during SWellEX-3, while the bottom portion of FIG. 14 shows therange-trajectory 1320 of the towed acoustic source. Note that the modelused reflects the initial portion of the bathymetry (up to time 60 min).After time 60 min., there is increasing mismatch between the model usedand the true propagation environment whose effect will be seen in thefinal portions of the tracks to be described next. An embodiment of themethods disclosed herein can alleviate this mismatch by dynamicallyupdating the acoustic propagation model used according to theappropriate bathymetric data along the direction on which the source islocated. This embodiment will use source location estimates (Eq. 3) todynamically modify the model-based predictor based on the current sourcelocations.

Referring to FIGS. 15-18, FIG. 15 shows a diagram 1400 illustratingdepth tracks obtained for one source using matched field tracking, FIG.16 shows a diagram 1500 illustrating depth tracks obtained for onesource using an embodiment of the methods disclosed herein, FIG. 17shows a diagram 1600 illustrating range tracks obtained for one sourceusing matched field tracking, and FIG. 18 shows a diagram 1700illustrating range tracks obtained for one source using an embodiment ofthe methods disclosed herein.

MFT used a square search window of 22 grid points in range and 5 gridpoints in depth. The MFT tracks are formed by the tracklets (shorttracks) that yielded the largest MFT score. The tracks corresponding toAlgorithm 2 show the range and depth of the 10-largest non-zero entriesin each SLM. Per SLM, the magnitude of all non-zero entries isnormalized with respect to the magnitude of the largest coefficientpresent. The number of non-zero rows of S, denoted S₀, obtained for eachSLM is shown above the tracks.

Despite its high computational complexity, MFT was used as a baselinefor constructing the source tracks. A total of 8 ambiguity surfacesobtained via Bartlett MFP, corresponding to 109 seconds of recordeddata, were used. Each ambiguity surface accounts for 13.65 seconds ofrecorded data. MFT tracks were incoherently averaged over frequency. ForAlgorithm 2, μwas adapted via grid search so as to obtain approximately10 nonzero entries per SLM. Thus, control on the sparsity of the trackswas exercised. Per time instant, the tracks were constructed by plottingthe range and depth of the largest coefficients in the correspondingSLMs. With our selection of tuning parameters, on average each SLM had 9nonzero entries. When all peaks were plotted to construct the tracks,some artifacts (horizontal lines) appeared in the tracks. Theseartifacts corresponded to source locations being maintained as part ofthe track and can be removed by using a dynamic selection scheme for λand μ.

By plotting fewer coefficients per SLM most of these artifacts can beremoved. Note that after time index 70, the source reaches an area withsignificantly different bathymetry. Thus, the difficulty that both MFTand the disclosed method have in tracking the source during the last legof the track is due to the mismatch between the environment and themodel used to construct the replicas. Note that a robust localizationframework can be used within the currently disclosed framework tomitigate the deleterious effect of model mismatch on the source-trackestimates. A specific embodiment of a robust localization framework thatcan be readily extended to the currently disclosed framework isdescribed in the commonly-assigned U.S. patent application Ser. No.14/285,400, entitled “Multitask Learning Method for BroadbandSource-Location Mapping of Acoustic Sources.”

FIG. 19 shows a diagram 1800 illustrating a Source Location Map (SLM)obtained using an accelerated proximal gradient (PG) solver. As shown,reference 1810 represents the estimated source location, which in thiscase matches the true source locations, while references 1820, 1830, and1840 represent artifacts, i.e., non-zero entries in the SLM that do notcorrespond to source locations. The color bar on the right illustratesthe intensity of each point on the map in a decibel (dB) scale. Notethat the magnitude of all artifacts is more than 10 dB below themagnitude of the point associated to the true source location.

FIG. 20 shows a graph 1900 illustrating faster convergence using anaccelerated PG solver versus a PG solver, where line 1910 represents thecost versus iterations for the accelerated PG solver and line 1920represents the cost versus iterations for the PG solver. In practice,only a few iterations of either algorithm were needed to correctlyidentify the support of the SLMs (up to 40 iterations in this case). Thedashed line 1930 illustrates the iteration value at which bothalgorithms were stopped, and the cost function value achieved by the PGand APG algorithms.

FIG. 21 shows a diagram 2100 illustrating range tracks obtained for twosources using an embodiment of a method in accordance with the PassiveTracking of Acoustic Sources with Sparse Innovations. In this case thepresence of two sources in the environment was created artificially bycombining two portions (first and second colored panels from left toright) of the one source track illustrated in FIG. 14. Similarly to theone source case, μwas adapted via grid search so as to obtainapproximately 15 nonzero entries per SLM. The top portion of FIG. 21illustrates the number of nonzero entries S₀ obtained per SLM. Thebottom portion of FIG. 21 illustrates the estimated trajectories for thetwo sources, which as expected match portions of the one sourcetrajectory shown in FIG. 14.

Although some embodiments of the method were discussed herein withregard to underwater source localization, some embodiments of the methodmay apply to other acoustic source localization environments, such asabove water, where accurate in-air acoustic propagation models areavailable. The embodiments of the method disclosed herein may also beextended to exploit prior information about source locations so as todevelop sparsity-cognizant tracking algorithms using passive sonar.Another possible extension involves using spatially-distributed arraysfor localization as a way to exploit spatial diversity to counteractmultipath affects in the localization performance and to reduce thepresence of surveillance gaps.

Many modifications and variations of the Passive Tracking of AcousticSources with Sparse Innovations are possible in light of the abovedescription. Within the scope of the appended claims, the embodiments ofthe systems described herein may be practiced otherwise than asspecifically described. The scope of the claims is not limited to theimplementations and the embodiments disclosed herein, but extends toother implementations and embodiments as may be contemplated by thosehaving ordinary skill in the art.

We claim:
 1. A method comprising the steps of: determining a grid

of G tentative locations r_(g) to localize more than one acousticsources; using an acoustic model to compute replicas that correspond tosimulated acoustic sources at locations in

, wherein the replicas comprise model-predicted Fourier coefficientvectors at F frequencies {ω_(f)}_(f=1) ^(F) corresponding to theacoustic pressure field as sampled by an acoustic sensor array having Nsensors; collecting, using the acoustic sensor array, time-series dataof actual acoustic measurements at each sensor of the acoustic sensorarray caused by the acoustic sources; using a short-time Fouriertransform (STFT) on the collected time series data partitioned into mblocks to compute Fourier coefficient estimates at frequencies{ω_(f)}_(f=1) ^(F) for all N sensors; constructing a set of Fouriercoefficient vectors Y^(m):=[y₁ ^(m), . . . , y_(F) ^(m)] ∈

^(N×F), where [Y^(m)]_(n,f) ∈

denotes the Fourier coefficient corresponding to ω_(f) for the n-thsensor in the m-th STFT block; modeling Fourier coefficient vectors atω_(f) for the m-th measurement block of the collected time series dataas y_(f) ^(m)=Σ_(g=1) ^(G)s_(g,f) ^(m)p_(g,f)+ε_(m,f), ∀f, where s_(g,f)^(m) denotes the unknown Fourier coefficient associated to the acousticsignature at frequency ω_(f) for a source located at r_(g), p_(g,f) ∈

^(N) is the replica for ω_(f) corresponding to a source located at r_(g)normalized so that ∥p_(g,f)∥₂=1, and ε_(m,f) denotes the Fouriercoefficients at ω_(f) corresponding to the noise in the m-th block;setting to zero all values of s_(g,f) ^(m) associated to tentativelocations r_(g)where acoustic sources are presumed to be absent;obtaining, at time τ_(m), an estimate S^(m) as the solution to${{\min_{S \in {\mathbb{C}}^{G \times F}}{\frac{1}{2}{\sum\limits_{f = 1}^{F}\; {{y_{F}^{m} - {P_{f}s_{f}}}}_{2}^{2}}}} + {\sum\limits_{g = 1}^{G}\; \left\lbrack {{\mu {ϛ_{g}}_{2}} + {\lambda {{ϛ_{g} - ϛ_{g}^{m - 1}}}_{2}}} \right\rbrack}},$where S:=[s₁, . . . , s_(F)] ∈

^(G×F),

′_(g) is the g-th row of S, P_(f):=[p_(1,f), . . . , p_(G,f)] ∈

^(N×G) is the matrix of replicas for ω_(f), and μ, λ>0 are tuningparameters; and generating one or more source location maps (SLMs) over

per frequency ω_(f) using S^(m), wherein each location on a particularSLM is associated with the magnitude of its corresponding acoustic gainestimate |ŝ_(g,f)], wherein estimates of the actual locations of the Kacoustic sources, {{circumflex over (r)}_(k)(τ_(m))}_(k∈κ), correspondto the locations of the K-largest coefficients |ŝ_(g,f)| depicted in theSLM.
 2. The method of claim 1, wherein the step of determining a grid oftentative locations of an acoustic source includes the step ofestimating the number of acoustic sources.
 3. The method of claim 1,wherein the step of generating SLMs over

per frequency ω_(f) comprises using each column s_(f) of S_(m) toconstruct the SLMs per frequency ω_(f).
 4. The method of claim 1,wherein the step of generating SLMs over

per frequency ω_(f) comprises generating a broadband SLM over

comprising all frequencies ω_(f) used to compute S^(m) using a whole rowof S^(m) for each r_(g) ∈

.
 5. The method of claim 1, wherein the step of generating SLMs over

per frequency ω_(f) comprises generating a broadband SLM over

using S^(m) by plotting the pairs (r_(g), ∥

_(g) ^(m)∥₂) for all r_(g) ∈

, where

_(g) ^(m):=[s_(g,1) ^(m), . . . , s_(g,F) ^(m)]′ ∈

^(F) comprises the entries of the g-th row of S^(m).
 6. The method ofclaim 1, wherein the estimates of the actual locations of the K acousticsources {{circumflex over (r)}_(k)(τ_(m))}_(k∈X) indexed by theindex-set X ⊂ {1, . . . , K} are given by the locations that correspondto the K-largest entries in the SLM, X ∈ max_(|x|=K) Σ_(κ∈X)∥

_(κ) ^(m)∥₂.
 7. The method of claim 1, wherein the SLMs are generatedusing an iterative problem solver based upon a proximal gradient method.8. The method of claim 1, wherein the acoustic sources are selected fromsurface-based acoustic sources and underwater acoustic sources.
 9. Themethod of claim 8, wherein the acoustic model is an underwater acousticmodel and the acoustic sensor array is a hydrophone array.
 10. A systemcomprising: a processor operatively connected to an acoustic sensorarray, wherein the processor is configured to perform the steps of:determining a grid

of G tentative locations r_(g) to localize more than one acousticsources; using an acoustic model to compute replicas that correspond tosimulated acoustic sources at locations in

, wherein the replicas comprise model-predicted Fourier coefficientvectors at F frequencies {ω_(f)}_(f=1) ^(F) corresponding to theacoustic pressure field as sampled by an acoustic sensor array having Nsensors; collecting, using the acoustic sensor array, time-series dataof actual acoustic measurements at each sensor of the acoustic sensorarray caused by the acoustic sources; using a short-time Fouriertransform (STFT) on the collected time series data partitioned into mblocks to compute Fourier coefficient estimates at frequencies{ω_(f)}_(f=1) ^(F) for all N sensors; constructing a set of Fouriercoefficient vectors Y^(m):=[y₁ ^(m), . . . , y_(F) ^(m)] ∈

^(N×F), where [Y^(m)]_(n,f) ∈

denotes the Fourier coefficient corresponding to ω_(f) for the n-thsensor in the m-th STFT block; modeling Fourier coefficient vectors atω_(f) for the m-th measurement block of the collected time series dataas y_(f) ^(m)=Σ_(g=1) ^(G)s_(g,f) ^(m)p_(g,f)+ε_(m,f), ∀f, where s_(g,f)^(m) denotes the unknown Fourier coefficient associated to the acousticsignature at frequency ω_(f) for a source located at r_(g), p_(g,f) ∈

^(N) is the replica for ω_(f) corresponding to a source located at r_(g)normalized so that ∥p_(g,f)∥₂=1, ε_(m,f) denotes the Fouriercoefficients at ω_(f) corresponding to the noise in the m-th block;setting to zero all values of s_(g,f) ^(m) associated to tentativelocations r_(g)where acoustic sources are presumed to be absent;obtaining, at time τ_(m), an estimate S^(m) as the solution to${{\min_{S \in {\mathbb{C}}^{G \times F}}{\frac{1}{2}{\sum\limits_{f = 1}^{F}\; {{y_{F}^{m} - {P_{f}s_{f}}}}_{2}^{2}}}} + {\sum\limits_{g = 1}^{G}\; \left\lbrack {{\mu {ϛ_{g}}_{2}} + {\lambda {{ϛ_{g} - ϛ_{g}^{m - 1}}}_{2}}} \right\rbrack}},$where S:=[s₁, . . . , s_(F)] ∈

^(G×F),

′_(g) is the g-th row of S, P_(f):=[p_(1,f), . . . , p_(G,f)] ∈

^(N×G) is the matrix of replicas for ω_(f), and μ, λ>0 are tuningparameters; and generating one or more source location maps (SLMs) over

per frequency ω_(f) using S^(m), wherein each location on a particularSLM is associated with the magnitude of its corresponding acoustic gainestimate |ŝ_(g,f)|, wherein estimates of the actual locations of the Kacoustic sources, {{circumflex over (r)}_(k)(τ_(m))}_(k∈κ), correspondto the locations of the K-largest coefficients |ŝ_(g,f)| depicted in theSLM.
 11. The system of claim 10, wherein the step of determining a gridof tentative locations of an acoustic source includes the step ofestimating the number of acoustic sources.
 12. The system of claim 10,wherein the step of generating SLMs over

per frequency ω_(f) comprises using each column s_(f) of S^(m) toconstruct the SLMs per frequency ω_(f).
 13. The system of claim 10,wherein the step of generating SLMs over

per frequency ω_(f) comprises generating a broadband SLM over

comprising all frequencies ω_(f) used to compute S^(m) using a whole rowof S^(m) for each r_(g) ∈

.
 14. The system of claim 10, wherein the step of generating SLMs over

per frequency ω_(f) comprises generating a broadband SLM over

using S^(m) by plotting the pairs (r_(g), ∥

_(g) ^(m)∥₂) for all r_(g) ∈

, where

_(g) ^(m):=[s_(g,1) ^(m), . . . , s_(g,F) ^(m)]′ ∈

^(F) comprises the entries of the g-th row of S^(m).
 15. The system ofclaim 10, wherein the estimates of the actual locations of the Kacoustic sources {{circumflex over (r)}_(k)(τ_(m))}_(k∈X) indexed by theindex-set X ⊂ {1, . . . , K} are given by the locations that correspondto the K-largest entries in the SLM, X ∈ max_(|X|=K)Σ_(κ∈X)∥

_(κ) ^(m)∥₂.
 16. The system of claim 10, wherein the SLMs are generatedusing an iterative problem solver based upon a proximal gradient method.17. The system of claim 10, wherein the acoustic sources aresurface-based acoustic sources.
 18. The system of claim 10 wherein theacoustic sources are underwater acoustic sources.
 19. The system ofclaim 10, wherein the acoustic model is an underwater acoustic model andthe acoustic sensor array is a hydrophone array.
 20. A systemcomprising: a processor operatively connected to a hydrophone arrayhaving N hydrophones, wherein the processor is configured to perform thesteps of: determining a grid

of G tentative locations r_(g) to localize more than one underwateracoustic sources; using an underwater acoustic model to compute replicasthat correspond to simulated acoustic sources at locations in

wherein the replicas comprise model-predicted Fourier coefficientvectors at F frequencies {ω_(f)}_(f=1) ^(F) corresponding to theacoustic pressure field as sampled by the N hydrophones; collecting,using the hydrophone array, time-series data of actual acousticmeasurements at each hydrophone of the hydrophone array caused by theunderwater acoustic sources; using a short-time Fourier transform (STFT)on the collected time series data partitioned into m blocks to computeFourier coefficient estimates at frequencies {ω_(f)}_(f=1) ^(F) for allN hydrophones; constructing a set of Fourier coefficient vectorsY^(m):=[y₁ ^(m), . . . , y_(F) ^(m)] ∈

^(N×F), where [Y^(m)]_(n,f) ∈

denotes the Fourier coefficient corresponding to ω_(f) for the n-thsensor in the m-th STFT block; modeling Fourier coefficient vectors atω_(f) for the m-th measurement block of the collected time series dataas y_(f) ^(m)=Σ_(g=1) ^(G)s_(g,f) ^(m)p_(g,f)+ε_(m,f), ∀f, where s_(g,f)^(m) denotes the unknown Fourier coefficient associated to the acousticsignature at frequency ω_(f) for a source located at r_(g), p_(g,f) ∈

^(N) is the replica for ω_(f) corresponding to an underwater acousticsource located at r_(g) normalized so that ∥p_(g,f)∥₂=1, and ε_(m,f)denotes the Fourier coefficients at ω_(f) corresponding to the noise inthe m-th block; setting to zero all values of s_(g,f) ^(m) associated totentative locations r_(g) where underwater acoustic sources are presumedto be absent; obtaining, at time τ_(m), an estimate S^(m) as thesolution to${{\min_{S \in {\mathbb{C}}^{G \times F}}{\frac{1}{2}{\sum\limits_{f = 1}^{F}\; {{y_{F}^{m} - {P_{f}s_{f}}}}_{2}^{2}}}} + {\sum\limits_{g = 1}^{G}\; \left\lbrack {{\mu {ϛ_{g}}_{2}} + {\lambda {{ϛ_{g} - ϛ_{g}^{m - 1}}}_{2}}} \right\rbrack}},$where S:=[s₁, . . . , s_(F)] ∈

^(G×F),

′_(g) is the g-th row of S, P_(f):=[p_(1,f), . . . , p_(G,f)] ∈

^(N×G) is the matrix of replicas for ω_(f), and μ, λ>0 are tuningparameters; and generating one or more source location maps (SLMs) over

per frequency ω_(f) using S^(m), wherein each location on a particularSLM is associated with the magnitude of its corresponding acoustic gainestimate |ŝ_(g,f)|, wherein estimates of the actual locations of the Kunderwater acoustic sources, {{circumflex over (r)}_(k)(τ_(m))}_(k∈κ),correspond to the locations of the K-largest coefficients |ŝ_(g,f)|depicted in the SLM.